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Abstract 



Gaussian approximations to the Boltzmann operator have proven themselves in recent years as useful tools for the study of the 
thermodynamic properties of rare gas clusters. They are, however, not necessarily correct at very low temperatures. In this article 
we introduce a first-order correction term to the frozen Gaus sian imaginary time propagator and apply it to the argon trimer. Our 
findings show that the correction term provides objective access to the quality of the propagator's results and clearly defines the 
"best" Gaussian width parameter. The strength of the correction monitored as a function of the temperature indicates that the 
1 5 results of the Gaussian propagator become questionable below a certain temperature. The interesting thermodynamic transition 
C*") from a bounded trimer to three body dissociation lies in the temperature range for which the Gaussian approximation is predicted 
to be accurate. 
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1. Introduction 

The imaginary time or Boltzmann operator exp( -f3H) be- 
longs to the most important quantities needed for understand- 
ing the equilibrium properties of multi-dimensional systems at 
finite temperatures. The derivatives of its trace give access to 
the mean energy and specific heat, which allow for the inves- 
tigation of thermodynamic properties of atomic clusters lUrO]. 
For systems with many degrees of freedom it is, however, still 
a challenge to evaluate the Boltzmann operator although it is 
accessible with Monte Carlo methods Il8l 4l0ll . since the neces- 
sary integrals can become numerically expensive, in particular 
at low temperatures. Approximations are essential to evaluate 
the Boltzmann operator for clusters of a few dozen atoms. 

Semiclassical initial value representations propagating Gaus- 
sian wave packets of the form 



<*|g> = (7r w |detG(r)|) 



•1/4 



x exp|--[* - q(T)] T G(TT l [x - q( T )] + y(r)|, (1) 



for the Boltzmann operator have been developed ll3l HI, 121 
and successfully applied to the study of thermodynamic prop- 
erties of clusters [3. In particular the time- 
evolved Gaussian approximation developed by Mandelshtam 
and co-workers has become an important tool which has 
been applied to thermodynamic properties of a large number of 
clusters J3I asm 11311 . This method is based on so-called 
thawed Gaussian wave packets, where the matrix of Gaussian 
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width parameters G(r) changes with (imaginary) time r. Re- 
cently a frozen Gaussian propagator with a constant width ma- 
trix G(r) = G(0) suggested by Zhang et al. lfl2Tl was shown 
to provide results competitive with the time-evolved Gaussian 
approximation using a single particle ansatz 01411 to simplify 
the computation of the time-dependent width matrix. Since no 
equations of motion for the N x N elements of G have to be 
propagated, the frozen Gaussian method is numerically much 
cheaper. 

Despite the large success and broad applicability of Gaus- 
sian methods for the evaluation of the mean energy and the spe- 
cific heat of clusters one has to keep in mind that they are semi- 
classical approximations. As pointed out by Liu and Miller [15] 
they are not expected to be accurate at very low temperatures 
where quantum effects are large. In the simplest case one ob- 
serves a shift of the ground state energy away from the exact re- 
sult [14J]. The low-temperature range is, however, very often the 
most important and crucial for the thermodynamic propertie s of 
rare gas cluster. For example, Gaussian based computations on 
clusters of light atoms, e.g., Nei3 an d Ne^s , predict novel 
low temperature quantum effects such as liquid-like zero tem- 
perature structures of Ne38 as compared to a solid-like structure 
predicted from classical mechanics IU6J . These predictions are 
based on the Gaussian approximation J3 ,|6 , 11 , 13 , 15 ] and have 
to be verified. 

A systematic approach connecting the Gaussian semiclassi- 
cal initial value approximations with exact quantum mechanics 
is the generalized time-depend ent p erturbation approach devel- 
oped by Pollak and coworkers 1 12, 17|. Within this framework 
the Gaussian approximations can be considered as a lowest- 
order approximation of a series converging to the exact quan- 
tum Boltzmann operator. It has been shown that both the thawed 
and frozen Gaussian versions of this series converge rapidly to 
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the numerically exact answer for a one-dimensional double well 
potential [12, 18]. In this paper we demonstrate that the correc- 
tions to the Gaussian imaginary time propagator are also prac- 
tically applicable to higher-dimensional systems, in particular 
to atomic clusters. To do so, we apply the frozen Gaussian se- 
ries to the partition function of the argon trimer and calculate 
its first-order correction. We show that the correction helps to 
estimate the quality of the results and provides objective access 
to the validity of the Gaussian approximations. With the first- 
order terms it is possible to identify a border temperature below 
which the Gaussian results become questionable. This result 
clearly states that the dissociation process ||7|] discussed earlier 
with Gaussian approximations [ 14] is correctly described by the 
frozen Gaussian imaginary time propagator since it appears in 
the allowed temperature range. 

The Gaussian width of the frozen Gaussian propagator is a 
free parameter but has an important impact on the accuracy of 
the physical quantities calculated with the approximation. It is 
demonstrated in this article that the first-order correction pro- 
vides an unambiguous method for optimizing the value which 
is based on minimizing the ratio between the first-order correc- 
tion to the zeroth-order term for the Boltzmann operator. This 
choice improves the results for the mean energy and the spe- 
cific heat. Including the first-order correction they reach the 
same quality as a fully-coupled thawed Gaussian computation 
over a wide range of temperatures. 

This article is organized as follows. In Sec.|2]we review the 
most important parts of the frozen Gaussian series represen- 
tation of the Boltzmann operator and introduce the first-order 
correction to the partition function for clusters of atoms. The 
theory is then applied to the argon trimer in Sec. [3] We intro- 
duce the system (Sec. 13. Il l, show how the first-order correction 
term may be used to determine the best width parameter (Sec. 
13.21 1. consider the influence of an artificial confinement on the 
correction (Sec. 13.3b and present the first-order corrected mean 
energy and specific heat for the dissociation process of the clus- 
ter (Sec. I3.4I >. which allow for a clear statement on the quality 
of the Gaussian approximation. The required numerical effort 
is analyzed in Sec.|4] Conclusions are drawn in Sec. [5] 

2. Frozen Gaussian approximation to the partition function 
and first-order corrections 

2.1. Frozen Gaussian series representation 

Zhang et al. Il2ll proposed a frozen Gaussian form of the 
imaginary time propagator K(J3) = exp(-/3H) at inverse tem- 
perature j6 = 1 / (kT) and developed its higher-order corrections 
in the framework of the generalized time-dependent perturba- 
tion series 1 17, T9L 20ll . The frozen Gaussian coherent state has 
the form 

(x\g(j,(r),q(T),r)) = {^^ 

X exp(-i[* - q(T)] T T[x - q(r)] + -V(r) • [x - 9 (t)A (2) 



with the symmetric positive definite matrix of constant width 
parameters T and the dynamical variables q{r) and p(t). Ap- 
proximating a solution to the Bloch equation 



—q-\Qo,t) = H\q ,r) 



(3) 



by propagating a frozen Gaussian wave packet, one finds that 
the imaginary time propagator matrix element has the form 



(x'\K (t)\x) = det(r)exp(--Tr(r)r 



^det(2 [1 - expt-^ror 1 ) 



xexp\--[x' - xYT[tanh(h £ TT/2)]- l [x' - x]\ 



'))> 



r da 3N i r /2 

x J (2^ exp rio dT ' mq{T,) 

- [x - <7(r/2)] T r[* - ? (t/2)] J (4) 

with x = (x' + x)/2. The only remaining dynamical variables 
in this representation of the propagator are the components of 
the vector q(j), which follow the equations of motion 



*£> = -r-<VV( 9 (r))>. 

OT 



(5) 



The angle brackets symbolize Gaussian averaged quantities of 
the form 

... .. /det(D\ 1/2 

/CO 
dx m exp(-[x-q] T r[x-q])h(x). (6) 

In the framework of the generalized time-dependent pertur- 
bation theory, Kq{t) is the zeroth-order approximation to the 
exact Boltzmann operator. Since the Bloch equation (0 is not 
solved exactly by the propagator (|4), one can define the correc- 
tion operator 01711 



C{t) = -—K (j)-HK {t). 

OT 



(7) 



The formal solution of Eq. (0 is 

K (t) = K(t)- f dT'K{T-T')C{r'). (8) 
Jo 

With the assumption that the exact propagator K(t) can be ex- 
pressed in terms of a series 



*(T) = 2*;(T), 



(9) 



n=0 



where Kj(t) ~ C(r) ; are terms with ascending power in the 
small correction, one can write down the recursion relation 



Kj +1 (T)= [ d T 'Kj(T-T')C(T') 

Jo 



(10) 



for the calculation of higher-order terms. 
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2.2. First-order corrections to the partition function 

In the framework of the series expansion the lowest-order 
approximation to the quantum partition function is simply the 
trace of the zeroth-order propagator ifTil . 

Z (y6) = Tr[£oG3)] = Vdet(r)exp(- jTrfl^J 



x A /det([l - exp(-h 2 T/3)Y l ) 
- dq m 



r°° dg 3N I r m \ 



Its first-order corrected counterpart 

Z 1 (/3)=Z (P)+Za(J3) 



(12a) 



is obtained in a similar way. The first-order correction to the 
propagator is evaluated according to Eq. (TTOt and then the cor- 
rection term for the partition function is 

Zb (0) = Tr [*!(£)] 
= J dT J dx' 3N j dx 3N (x'\K (J3-T)\x)(x\C(T)\x'). (12b) 

More details of the correction operator and the first-order term 



to the partition function can be found in Appendix A In par- 
ticular it is shown, how the x' and x integrations in Eq. (I12bb 
can be performed analytically if the potential of the system can 
be expressed in terms of Gaussians. This makes the evaluation 
of the correction terms efficient. 

Note that the first-order correction term (112b) is not sym- 
metric. Here we used the "left" correction operator defined by 
Zhang et al. Il2ll . This has, however, no influence on the parti- 
tion function due to the trace in Eq. (112bl i. 

3. First-order calculations for the argon trimer 

3.1. Potential and confinement 

To be consistent with previous studies of the system 10, EH 
2111 we model the pairwise interaction of the argon atoms by a 
Morse potential 

V(r u ) = D (exp [-2a(r u - R e )] - 2 exp [-o(ry - R e )]) ( 1 3) 

with the distance r (J between particles i and j and the Morse 
parameters D = 99.00cm" 1 , a = 1.717 A, and R e = 3.757 A. 
Since the Gaussian averaged quantities according to Eq. © and 
the x' and x integrations in the first-order correction term (112bt 
can be performed analytically if the potential can be expressed 
in terms of Gaussians, we fitted the Morse potential (TT3l to a 
sum of three Gaussians 

3 

V(|r, - rj\) = c p e- a ^, r u = \n - r,|, (14) 
with the parameters listed in Table Q] 



Table 1 : Parameters used in the Gaussian representation d!4t of the potential, 
taken from Ret'. | 



p 


c p [cm '] 


a p [A 2 ] 


1 


3.296 x 10 5 


0.6551 


2 


-1.279 x 10 3 


0.1616 


3 


-9.946 x 10 3 


6.0600 



Since all x' and x integrations in the zeroth-order approxi- 
mation and the first-order correction term of the propagator can 
be performed analytically with the Gaussian potential (TT41 only 
the integration in the q space and the imaginary time integra- 
tion in Eq. d!2bl ) remain for numerical computation. The t in- 
tegration is one-dimensional and can be done with a standard 
integrator. The multi-dimensional Monte Carlo sampling in the 
q variables is done with a standard Metropolis algorithm as out- 
lined in Ref. H. 

To converge the integrals the q space is restricted by the 
condition \q - R cm \ < R c , where R c is a confining radius and 
R cm is the center of mass of the cluster. That is, no atom may 
leave the center of mass beyond a certain distance R c . It is also 
possible to introduce the confinement via an additional steep 
potential as 



V c (r) * 



Rc 



20 



(15) 



As has been shown previously , the value of R c has 

a critical influence on the thermodynamic properties. For the 
first-order correction term (112b) two q integrations, viz. in the 
propagator and in the correction operator, have to be taken into 
account. Here, the restriction on q is applied for the propagator 
part (@). Due to the coupling via x' and x the restriction auto- 
matically affects also the q variable in the correction operator. 



More details can be found in Appendix A 



In this article we concentrate on a confinement of R c = 
10 A. It is known that this confinement does not describe the 
dissociation process fully since it is too restrictive. However, 
the restricted cluster shows a very broad transition from the 
bounded system to three free particles, and thus allows for a 
clear analysis of the correction term's influence on the results 
around the dissociation, which is the interesting physical effect. 

3.2. Choice of the width parameter 

In the frozen Gaussian ansatz the matrix of width parame- 
ters r is constant in imaginary time t, however, the choice of 
r has a critical influence on the quality of the results lfl2ll . In a 
previous study [14] we found that the width matrix 



r = 



((D 1+ 2D 2 )/3 
(Di-D 2 )/3 
(Di-D2)/3 



(d - D 2 )/3 
(D l+ 2D 2 )/3 
(Di-Da)/3 



(D l -D 2 )/3 
(D 1 - D 2 )/3 
(D 1+ 2D 2 )/3) 



(16) 



with two 3 x 3-submatrices D\ — D\ \ and D 2 = D 2 1, where 
1 is the 3x3 identity matrix, provides results competitive with 
thawed Gaussian calculations using a single-particle ansatz. The 
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parameter D\ represents the free center of mass motion and 
should be chosen as small as possible. The Gaussian width for 
the relative coordinates is represented by Do- In Ref. IU4I1 an 

o -2 

optimum value of Dx — 25 A was found. It yielded the best 
approximation to a numerically exact path integral Monte Carlo 
calculation [7] and thawed Gaussian approximations lfl4ll . In 
the same work it was found that all choices of D\ below a cer- 

• -2 

tain value lead to the same mean energies and that D\ — 0. 1 A 
belongs to this range in which D\ is small enough. These pa- 
rameters also provide the best and lowest approximation to the 
ground state energy. Thus, the minimum of the ground state en- 
ergy could also be used to search for the optimum value of the 
width parameters. It provides, however, no access to the quality 
of the results if the exact ground state energy is unknown. 

A systematic, objective and internally consistent approach 
to the determination of the optimum width parameter is made 
possible by considering the first-order correction. The defini- 
tion of the correction operator in Eq. (|7} makes clear that it 
measures the deviation of the zeroth-order term from the exact 
result. One should thus minimize the ratio of the first-order [Eq. 
d!2bH and zeroth-order [Eq. (ITTt l contributions to the partition 
function by varying the width parameter matrix and thus find 
the optimum value. 

In Fig. Q2a) we plot the ratio Zqi/Zq for several values of the 

° -2 

parameters D2 and D\ — 0.1 A . The data was calculated for a 
cluster with a confinement of R c = 10 A. The zeroth- and first- 
order mean energies in Fig. [TJb) visualize the transition from 
a bounded to an unbounded cluster. For temperatures at which 
the system is bounded, i.e., T ig 30 K an optimum range for the 
width parameter at which the contribution of the correction term 

Zci is minimal can be found. Below the transition temperature 

0-2 0-2 

the lines for D2 between 20 A and 30 A are close to each 

-2 

other, the smallest value is found for D2 — 25 A . This con- 



firms our previous zeroth-order study of the argon trimer [14], 

° -2 

where the same value of D2 — 25 A provided the best agree- 
ment with numerically exact path integral Monte Carlo meth- 
ods. At higher temperatures, where the cluster passes through 
a transition to three almost free atoms, the situation changes. 
The ratio becomes smaller as the magnitude of D2 decreases. 
In the limit of three free particles, the frozen Gaussian approxi- 
mation is exact, in the limit that the width matrix vanishes. One 
thus expects that in this limit, the smaller values of the width 
parameters would provide a better approximation. 



3.3. Influence of the confining sphere 

As already mentioned, for free particles, the frozen Gaus- 
sian propagator yields the exact partition function in the limit 
r — > 0. Thus, one expects that the free center of mass motion, 
described by D\, demands D\ to be as small as possible. In 
practice one observes, however, that D\ shows a minimum cor- 
rection at a finite value. This is a result of the confinement of the 
atoms, which has the same effect as adding a confining poten- 
tial ( fTBI l, i.e., the particles are not really free. The effect can be 
observed in Fig. [2] where the ratio Zci /Zo is plotted for several 

o -2 

choices of D\ and D2 — 25 A .As can be seen the minimum 
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Figure 1 : (a) The ratio of the contribution of the first- and zeroth-order terms 
in the series expansion of the partition function is plotted as a function of the 
temperature T for varying values of the width parameter D2 (given in units of 

A~ ) and D\ = 0.1 A" . For the bounded system (T < 30 K) the range D 2 = 

° -2 ~ 
20 ... 30 A leads to the minimal ratio whereas for temperatures above the 

transition temperature smaller values of Di lead to smaller relative corrections, 
(b) The zeroth-order approximation to the mean energy Eq and its first-order 
corrected counterpart E\ is shown to visualize the transition from the bounded 

9-2 j -2 

to the unbounded system for D\ = 0. 1 A and D2 = 25 A . 

" -2 

correction is achieved for the finite values D\ = 0.5 ... 1 A if 
a confinement of R c = 10 A is used. 

The influence of the confinement can also be found in the 
study of the width parameter D2 for the internal degrees of free- 
dom. For high temperatures, i.e., in the limit of three free par- 
ticles a smaller value of D2 should provide a better approxima- 
tion. Thus, one expects always a smaller contribution Zci for 
smaller values of D2. However, the strengths of the first-order 

0-2 o-2 

corrections for Z>2 = 10 A and D2 — 15 A are almost the 
same in the high-temperature limit, even beyond the tempera- 
tures shown in the figure. 

3.4. Correction at low temperatures and the first-order cor- 
rected mean energy 

Figures 02a) and [2] already indicate that the correction term 
increases in importance as the temperature is lowered. This 
is not surprising since an imaginary time propagation is per- 
formed approximately and the difference between the approx- 
imate and exact solutions is expected to increase with (imagi- 
nary) time. The behavior becomes even clearer in Fig. [3] where 
the relative strength of the first-order correction Zc\/Zq is plot- 
ted for the argon trimer with a confinement R c = 10 A and width 
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Figure 2: The ratio of the contribution of the first- and zeroth-order terms in 
the series expansion of the partition function is plotted as a function of the 
temperature T for varying values of the width parameter D\ (given in units 

9 -2 9 -2 

of A ) and Di = 25 A . The minimal ratio is found in the range D\ = 

a -2 

0.5 ... 1 A . The cluster is limited by a confinement with R c = 10 A. 
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Figure 3: Relative strength of the first-order correction to the partition function 

o-2 9-2 
at lower temperatures for R c = 10 A, D\ = 0.1 A , and D2 = 25 A .It 

diverges in the limit T — > 0. 



o-2 0-2 

parameters Di = 0.1 A , D2 = 25 A down to a temperature 
of T — IK. It increases drastically and shows a diverging be- 
havior in the limit T — » 0, i.e., /3 — > 00. The increasing impor- 
tance of the first-order correction term at low temperatures is 
consistent with our previous zeroth-order study of the system, 
where we found that the accuracy of the frozen Gaussian ap- 
proximation gets worse as T —> and the cluster approaches 
the ground state II 1 411 . At low temperatures, quantum effects are 
strong and, in particular, the Gaussian (zeroth-order) approx- 
imations cannot reproduce the exact ground state energy £gs 
of the cluster. The Gaussian form imposed on the wave func- 
tion is too severe and the ground state energy obtained from the 
zeroth-order approximation is found to be slightly higher than 
the exact ground state energy (£0 = £gs + AZs and AE > 0). 

To further understand the divergence of the first-order term, 
we may assume that for large /3 only the ground state contributes 
to the partition function and that the exact partition function Z 
and its zeroth-order approximated counterpart Zq behave as 



Z 



-/3(E as +AE) 



(17) 



■frozen Gaussian, zeroth order — 
frozen Gaussian, first order - • 
thawed Gaussian 




Figure 4: (a) Temperature dependence of the mean energy of the argon trimer 

9 -2 9-2 
for R c = 10 A, D| = 0.1 A , and D2 = 25 A . Shown are the zeroth-order 



(solid line) and first-order (dashed line) frozen Gaussian results, and the fully- 
coupled thawed Gaussian approximation (dotted line), (b) Specific heat for the 
same parameters. Both calculations indicate a substantial improvement due to 
the first-order correction in the temperature region of the transition and above. 
For lower temperatures the correction breaks down, indicating that the frozen 
Gaussian approximation is no longer sufficiently accurate. 



Thus, a first-order correction which is valid also for large j3 
would be expected to diverge, since Z/Zo ~ e^ AE . 

More interesting is, however, the influence of the first-order 
correction on the derivatives of the partition function, which are 
required for the physically meaningful average energy 



E = kT 2 



dlnZ 
dT 



and specific heat 



dT 



(18a) 



(18b) 



Figure [4] shows the temperature dependence of the mean en- 
ergy and specific heat of the argon trimer for R c = 10 A, D\ = 

9 — 2 9 — 2 

0.1 A , and D2 — 25 A , i.e., the same parameters as in Fig. [3] 

The frozen Gaussian results for the zeroth- and first-order ap- 
proximations are compared with a fully-coupled thawed Gaus- 
sian approximation which, as discussed elsewhere lfl4l is ex- 
pected to be more accurate then the frozen Gaussian approxi- 
mation. As one can see in Fig. |4fa) the first-order correction 
for the frozen Gaussian approximation brings the resulting esti- 
mate closer to the fully-coupled thawed Gaussian estimate over 
most of the temperature range studied. However, for tempera- 
tures T ;$ 12 K the first-order correction to the energy becomes 
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Figure 5: Double logarithmic plot of the difference D between the first deriva- 
tives of the logarithms of the zeroth- and first-order partition functions as de- 
fined in Eq. )19t . In the limit of low temperatures the divergence behaves as 
D ~ 7"~ - 91 , which is not strong enough to compensate for the T 2 term in Eq. 
(18a) needed to obtain a nonvanishing energy correction. 



smaller and vanishes in the limit that T — » 0. This also in- 
fluences the specific heat as can be seen in Fig. Sfb). In the 
temperature range of the transition from bounded to unbounded 
motion, the first-order corrected frozen Gaussian calculations 
agree well with the fully-coupled thawed Gaussian counterpart. 
For lower temperatures this is no longer so. The divergence is 
especially noticeable in the low temperature limit for the spe- 
cific heat where the first-order correction significantly increases 
the difference between the fully-coupled thawed Gaussian esti- 
mate (which is rather accurate as known from numerically exact 
computations, see Ref. UH) and the frozen Gaussian based es- 
timate. 

Interestingly, although the first-order correction term at very 
low temperatures is much larger than the zeroth-order frozen 
Gaussian estimate (see Fig. [3]), it does not lead to any change in 
the estimate of the ground state energy. As may be seen from 
Fig- HI a) the increase of the correction to the partition function 
is not "fast enough". This is further demonstrated in Fig. [5] 
where the difference between the first derivatives of the loga- 
rithms of the zeroth-order and first-order partition functions, 



D = 



d In Zq d In Z\ 



3T 



8T 



(19) 



is plotted vs. the temperature. From Fig.[3]one expects this dif- 
ference to diverge and this is evidently the case. However, the 
first-order correction term can lead to a correction for the e sti- 
mate of the ground state energy only if the temperature deriva- 
tive diverges as T 2 , as can be seen from Eq. (II 8ab . We find, 
however, that below a critical temperature the difference of the 
derivatives behaves as D ~ T~ 91 , and this is not sufficient. 
As a result the first-order correction to the ground state energy 
vanishes as found in Fig. [4] 

Although the first-order correction term does not lead to a 
corrected estimate for the ground state energy, it does provide 
information on the quality of the frozen Gaussian approxima- 
tion. The ratio Zci/Zq (cf. Fig. [3j> and the temperature behav- 
ior of the correction to the mean energy [cf. Figs. 0(a) and|5] 
provide an objective measure for the correctness of the values 



as estimated from the frozen Gaussian approximation. Only 
at a temperature T ~ 12 K and below it, the correction be- 
comes comparable with the zeroth-order estimate, and thus is 
no longer small. Approximately at the same temperature the 
energy correction starts to get smaller. Both outcomes indicate 
that the quality of the approximation is questionable for lower 
temperatures. However, the dissociation process described in 
Ref. [14] appears only at higher temperatures, indicating that 
the frozen Gaussian approximation accurately reflects the quan- 
tum transition from a bounded moiety to dissociation into three 
free particles. 

4. Numerical effort for evaluating the first-order correction 

The results presented for the mean energy and the specific 
heat in Sec. [3]showed that the quality of the first-order corrected 
frozen Gaussian approximation is similar to the fully-coupled 
thawed Gaussian method in the interesting range of tempera- 
tures. It is therefore of interest to compare their computational 
cost. 

The numerically most expensive part is the imaginary time 
propagation of the dynamical variables. Since the positions q 
in the zeroth-order frozen Gaussian propagator (3] and the cor- 
rection operator ([7} are the only dynamical variables of a frozen 
Gaussian propagator, the number of equations of motion for the 
first-order corrected partition function still scales linearly with 
the particle number. For N particles one has to solve 6N + 2 
equations of motion, in which the time integrations of the aver- 
aged potential in the exponential of the propagator and the cor- 
rection operator [cf. Eq. (|4} and Appendix A] are included. By 
contrast the number of dynamical variables of the fully-coupled 
thawed Gaussian propagator scales always quadratically with 
the number of particles due to the symmetric time-dependent 
Gaussian width matrix. Counting all dynamical variables of 
the time-evolved Gaussian approximation J3, IT] one obtains 
3N(3N + 3)/2 + 1 equations of motion. Thus, the imaginary 
time propagation of the dynamics is cheaper for the first-order 
frozen Gaussian partition function. 

The first-order correction requires, however, also additional 
numerical effort. The need for propagating pairs of q trajecto- 
ries for the propagator and the correction operator in Eq. (II 2b\ 
demands a larger number of sampling points for the Monte 
Carlo integration in the two sets of positions q. Additionally, 
the time integration in Eq. ( 1 12bb has to be evaluated. At this 
point the importance of an analytical evaluation of the x' and x 
integrations in Eq. (I12bb as described in Appendix A becomes 
clear. Owing to this simplification one can do without an ad- 
ditional 6Af-dimensional Monte Carlo sampling for the opera- 
tor product [x integration in Eq. (112bH and the trace (x' inte- 
gration). This drastically reduces the numerical costs. Note 
that a 3Af-dimensional trace integration is already required for 
the zeroth-order approximation of both the frozen and thawed 
Gaussian partition functions. It can be and is evaluated ana- 
lytically in all cases. With the analytical evaluation of the x' 
and x integrations there remain 6N q integrations for the first- 
order corrected frozen Gaussian partition function and 3iV for 
its zeroth-order thawed Gaussian counterpart. 
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Because of the different numbers of required sampling points 
there is no general way for estimating the numerical effort. In 
the case of the argon trimer the additional effort for the first- 
order correction outweighs the lower number of equations of 
motion. We used 6.5 x 10 7 sampling points for the frozen Gaus- 
sian with first-order correction and 1.2 x 10 7 for the thawed 
Gaussian to obtain converged results down to low temperatures 
of about 1 K. On the same architecture (one Tesla CI 060 GPU) 
the frozen Gaussian first-order calculation resulted in approx- 
imately twice the time as the zeroth-order thawed Gaussian. 
Due to the better scaling of the number of dynamical variables 
(linear instead of quadratic) one may expect that this relation 
changes in favor of the first-order frozen Gaussian variant with 
increasing particle number. However, as the results in this arti- 
cle demonstrate, the information the first-order correction pro- 
vides on the validity of the Gaussian approximations is at least 
as important as the improvement of the numerical estimate of 
the physical values. 



confirmed with the correction term and whether the breakdown 
of the approximation lies at a lower temperature than the frozen 
Gaussian results. 

Since most effects in rare gas clusters such as structural 
transformations or dissociations appear at low temperatures fl- 
it will be of value to investigate these properties with the 
series expansion. The correction term should be used to verify 
the validity of the Gaussian approximations in these cases, in 
particular, where strong differences are found between the ap- 
proximate quantum computations and a purely classical theory 

m 
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5. Conclusions and outlook 

In this article we introduced first-order corrections to a frozen 



Gaussian approximation [12] of the Boltzmann operator with 
application to the thermodynamics of atomic clusters. By us- 
ing a Gaussian fit of the underlying potential it is possible to 
evaluate many of the necessary integrations of the correction 
terms analytically such that the first-order correction becomes 
viable for systems with "many" degrees of freedom. We ap- 
plied the correction to a study of the thermodynamics of the 
argon trimer, whose dissociation process has been investigated 
recently iBB- 

The highest value of the correction term is to asses the qual- 
ity of Gaussian approximations used in the study of thermody- 
namic properties of high-dimensional systems llll-l7[[l3, 14], It 
is known that these Gaussian approximations are exact in the 
high-temperature limit but are not necessarily correct at low 
temperatures. The present study shows that the first-order cor- 
rection indicates a border temperature below which the results 
of the approximate propagators become questionable and above 
which they may be considered to be reliable. 

The investigation of the argon trimer revealed that the dis- 
sociation process J3l discussed earlier with Gaussian approxi- 
mations in full detail II 1 411 is correctly described by the frozen 
Gaussian imaginary time propagator. It appears in the temper- 
ature range for which the first-order correction is small. Fur- 
thermore, the first-order corrected results are comparable with 
a fully-coupled thawed Gaussian investigation of the system, 
i.e., the addition of the first-order correction term improves the 
thermodynamic estimates. The numerical cost of the first-order 
corrected frozen Gaussian values is, however, higher than that 
of the thawed Gaussian partition function. 

The series expansion of the Boltzmann operator exists also 
for thawed Gaussian propagators [17] and has been applied to a 
one-dimensional system [18]. Since the thawed Gaussian prop- 
agator approximates the exact mean energy of the argon trimer 
at low temperatures better than its frozen Gaussian counterpart 
discussed here, it will be interesting to see whether this can be 



Appendix A. First-order correction term to the partition 
function for a potential expressed in terms of 
Gaussians 

The correction operator for the frozen Gaussian imaginary 
time propagator can be written as [12] 



<*'|C(r)|*> = <*'| [~ - HjK (T)\x) = det(D 

X exp|-^-Tr(r)-rj ^det(2 [1 - exp(-fi 2 rr)] _I ) 
x expl ~[x' - *] T r[tanh(ft 2 rT/2)rV - x]\ 



r da 3N i r' 2 

x J ^—AV(x',x,q(T/2))ex V \-2 ^ dr' {V{q(r'))) 

-[jc- q(T/2)] T T[x - q(T/2)] J, (A. 1 a) 
where the energy difference operator is found to be 

AV(x',x, q(r/2)) = ^ J-Tr(D + [*' - *] T y [*' - x] J 

+ <V(?(T/2))> + j([x - q(T/2)] J T 2 [x - q(r/2)] 

+ [*' -*] T rcom(^rJr[jc - ? (r/2)]J 

+ [x - q(r/2)] ■ (VV(q(r/2))) - V(x'). (A. lb) 

Since most parts of the propagator <j4j and of the correction 
operator (lA.lat include only simple Gaussians or polynomials 
in x and x' the matrix element of the product of the zeroth-order 
propagator and the correction operator needed for computation 
of the first-order correction term may be calculated analytically. 
Only the parts including the potential require in general a nu- 
merical computation. However, assuming a Gaussian form (fT4b 
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for the potential enables one to calculate this part analytically 
as well. That is, all x' and x integrations can be performed 
analytically and only the q and r integrations remain for the 
numerical evaluation. After evaluating the x' and x integrations 
in Eq. d!2bb we can write the first-order correction as 

Zci(/?)= i2^ exp r- Tr(r) ^j 

x dr A /det([l - exp(-h 2 T(J3 - t))\ x ) 

x ^det([l - exp(-?i 2 rT)] _1 ) 

x j dq\ N j dq\ N AV s {q u q 2 ,^T) 

I H/3-t)/2 

xexpl-2 J o dT'(V( qi {T'))) 
xexp(-2 dT'<F( 92 (r')))), 



(A.2) 



where the single trajectory contribution to the energy difference 
operator is denoted as 



A\' S ( </ (li.fi. t) = exp |-~A?(/3, T) T TAq(/3, r) 



Vdet(B(y3,r))det(2r) 



+ ^Tr(B(J3, r)-'r 2 ) + l -Aq(J3, T) T T 2 Aq(JJ, r) 
+ (V(q 2 (T/2))) + l -Aq{p,T) T {VV{q 2 {rl2))) 



-Vsa{q,p,T)\, (A.3) 



the difference trajectory is defined as 

Aq(B, t) = f i(G8 - r)/2) - q 2 (r/2), (A.4a) 
and the middle point is 

m r) = \{qW ~ t)/2) + ? 2 (r/2)). (A.4b) 

The integrated potential Vsg(?,jS, t - ) can be written in a form 
similar to the Gaussian average of the Gaussian fitted potential 
( fT4b . which was introduced by Frantsuzov et al. ]3|] for a fully- 
correlated multi-dimensional Gaussian propagator. It consists 
of a sum of Gaussians, 



Vso(§,j8,t): 



det(CQB,T)) 
^det(B( y 6,T) + I72) 



EE' 



det(fl, 7 0g,T)) 
^det(Dy08,T) + ar J ,) 



: exp (- [^08, r) - ^-03, r)] T Fj^, t) r) - ^(JS, r)]) , 

(A.5) 



with the matrices 

B{fi, t) = ir [tanh(ft 2 r(j6 - t)/2) _1 + tanh(ft 2 rr/2)- 1 ] 

(A.6a) 



and 



C(J3,T) = UB(J3,T) + ^)[B(J3,T)Tr 1 



(A.6b) 



The submatrices C,- 7 -0S, r) for the particles i and y of C(J3, r) are 
required for 

D u (fi, t) = (CsOS, r) + Cjj(J3, t) - Cy(j8, r) - C^, r))"' 

(A.7a) 



and 



FgCS, t) = a p - a 2 p (a p + D U (J3, t)\ 



(A.7b) 



In Eq. (I A. 31 > one can clearly see that due to the Gaussian 
in A^ the contribution of a pair of trajectories q i(t) and q%(j) 
will only be nonvanishing if the distance between them is small. 
This implies that the restriction of q\ to the confining sphere 
defined by R c will automatically impose a restriction on q 2 . 
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